Preprocessing

library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)

df <- read.csv("../data/study1.csv") |>
  mutate(
    Date = as.Date(Date),
    Participant = fct_reorder(Participant, Date),
    Screen_Refresh = as.character(Screen_Refresh),
    Illusion_Side = as.factor(Illusion_Side),
    Block = as.factor(Block),
    Education = fct_relevel(Education, "Master", "Bachelor", "High School", "Other")
  )

Exclusions

outliers <- c(
  # Half of the trials have very short RT
  # Prolific Status: REJECTED
  "S33",
  # Block n2 with very short RTs
  # Prolific Status: RETURNED
  "S20",
  # Error rate of 46.2% and short RTs
  # Prolific Status: RETURNED
  "S51",
  # Error rate of 46.2%
  # Prolific Status: RETURNED
  "S49",
  # Error rate of 47.9%
  # Prolific Status: REJECTED
  "S47",
  # Error rate of 42.1% and very large RT SD
  # Prolific Status: REJECTED
  "S12"
)

We removed 6 participants upon inspection of the average error rage (when close to 50%, suggesting random answers) and/or when the reaction time distribution was implausibly fast.

For each block, we computed the error rate and, if more than 50%, we discarded the whole block (as it likely indicates that instructions got mixed up, for instance participants were selecting the smaller instead of the bigger circle).

Error Rate

dfsub <- df |>
  group_by(Participant) |>
  summarize(
    # n = n(),
    Error = sum(Error) / n(),
    RT_Mean = mean(RT),
    RT_SD = sd(RT),
  ) |>
  ungroup() |>
  arrange(desc(Error))

knitr::kable(dfsub) |> 
  kableExtra::row_spec(which(dfsub$Participant %in% outliers), background  = "#EF9A9A")
Participant Error RT_Mean RT_SD
S47 0.479 262 296
S51 0.465 263 372
S49 0.462 630 611
S12 0.421 507 1679
S20 0.402 492 725
S41 0.300 723 579
S50 0.287 659 333
S16 0.269 578 172
S05 0.262 716 271
S04 0.262 588 206
S42 0.260 718 1294
S43 0.255 837 711
S33 0.251 599 1326
S46 0.246 699 414
S06 0.246 560 197
S39 0.243 693 232
S15 0.238 1160 1660
S40 0.227 748 243
S08 0.225 780 427
S27 0.219 785 836
S30 0.216 506 150
S26 0.215 738 375
S10 0.213 1076 557
S23 0.210 870 489
S52 0.207 804 378
S34 0.206 652 367
S28 0.205 511 147
S32 0.203 627 223
S35 0.202 1133 982
S44 0.202 869 424
S24 0.201 720 298
S01 0.198 1110 738
S48 0.196 867 652
S18 0.189 983 830
S45 0.187 711 195
S25 0.186 803 397
S36 0.185 846 765
S09 0.182 1045 1082
S19 0.181 1001 562
S22 0.176 1307 1091
S14 0.175 1146 900
S02 0.174 703 243
S29 0.173 1046 918
S38 0.173 616 381
S17 0.171 927 550
S37 0.170 769 312
S13 0.157 848 364
S31 0.156 712 383
S07 0.142 1109 863
S11 0.139 895 816
S03 0.124 741 331
S21 0.092 884 509

Error Rate per Illusion Block

temp <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  summarize(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |> 
  arrange(desc(ErrorRate_per_block))

temp2 <- temp |> 
  filter(ErrorRate_per_block >= 0.5) |> 
  group_by(Illusion_Type, Block) |> 
  summarize(n = n()) |> 
  arrange(desc(n), Illusion_Type) |> 
  ungroup() |> 
  mutate(n_trials = cumsum(n * 56),
         p_trials = n_trials / nrow(df))

# knitr::kable(temp2)

p1 <- temp |>
  estimate_density(at = c("Illusion_Type", "Block")) |>
  ggplot(aes(x = x, y = y)) +
  geom_line(aes(color = Illusion_Type, linetype = Block)) + 
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(y = "Distribution", x = "Error Rate") +
  theme_modern()

p2 <- temp2 |> 
  mutate(Block = fct_rev(Block)) |> 
  ggplot(aes(x = Illusion_Type, y = p_trials)) +
  geom_bar(stat="identity", aes(fill = Block)) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
  labs(y = "Percentage of Trials Removed", x = "Illusion Type") +
  theme_modern() +
  theme(axis.text.x = element_text(angle=45, hjust = 1))

p1 | p2



# Drop
df <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  mutate(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |> 
  filter(ErrorRate_per_block < 0.5) |>
  select(-ErrorRate_per_block)

rm(temp, temp2)

Reaction Time Distribution

# RT distribution
p <- estimate_density(df, select = "RT", at = c("Participant", "Block")) |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  mutate(color = ifelse(Participant %in% outliers, "red", "blue")) |>
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = color, group = interaction(Participant, Block), linetype = Block)) +
  geom_vline(xintercept = 2500, linetype = "dashed", color = "red") +
  scale_color_manual(values=c("red"="red", "blue"="blue"), guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(xlim = c(0, 3000)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  facet_wrap(~Participant) +
  labs(y = "", x = "Reaction Time (ms)")
p

# ggsave("figures/outliers.png", p, width=20, height=15)

# Filter out
df <- filter(df, !Participant %in% outliers)

Reaction Time per Trial

p1 <- estimate_density(df, select = "RT", at = "Participant") |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = Participant, group = Participant)) +
  geom_vline(xintercept = c(150, 3000), linetype = "dashed", color = "red") +
  scale_color_material_d("rainbow", guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(xlim = c(0, 3500)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  # facet_wrap(~Participant) +
  labs(y = "", x = "Reaction Time (ms)")


df$Outlier <- df$RT < 150 | df$RT > 3000

p2 <- df |>
  group_by(Participant) |>
  summarize(Outlier = sum(Outlier) / n()) |>
  mutate(Participant = fct_reorder(Participant, Outlier)) |>
  ggplot(aes(x = Participant, y = Outlier)) +
  geom_bar(stat = "identity", aes(fill = Participant)) +
  scale_fill_material_d("rainbow", guide = "none") +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
  see::theme_modern() +
  theme(axis.text.x = element_blank())

p1 | p2

We removed 692 (1.37%) outlier trials (150 ms < RT < 3000 ms).

df <- filter(df, Outlier == FALSE)

Participants

dfsub <- df |>
  group_by(Participant) |>
  select(Participant, Age, Sex, Education, Nationality, Ethnicity, Duration, Break_Duration, Screen_Resolution, Screen_Refresh, Device_OS) |>
  slice(1) |>
  ungroup()

The final sample included 46 participants (Mean age = 26.7, SD = 7.7, range: [19, 60]; Sex: 39.1% females, 56.5% males, 4.3% other; Education: Master, 17.39%; Bachelor, 41.30%; High School, 39.13%; Other, 2.17%).

plot_distribution <- function(dfsub, what = "Age", title = what, subtitle = "", fill = "orange") {
  dfsub |>
    ggplot(aes_string(x = what)) +
    geom_density(fill = fill) +
    geom_vline(xintercept = mean(dfsub[[what]]), color = "red", linetype = "dashed") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    ggtitle(title, subtitle = subtitle) +
    theme_modern() +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      plot.subtitle = element_text(face = "italic", hjust = 0.5),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank()
    )
}

plot_waffle <- function(dfsub, what = "Nationality") {
  ggwaffle::waffle_iron(dfsub, what) |>
    # mutate(label = emojifont::fontawesome('fa-twitter')) |>
    ggplot(aes(x, y, fill = group)) +
    ggwaffle::geom_waffle() +
    # geom_point() +
    # geom_text(aes(label=label), family='fontawesome-webfont', size=4) +
    coord_equal() +
    ggtitle(what) +
    labs(fill = "") +
    theme_void() +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
p1 <- plot_distribution(dfsub, "Age", fill = "#FF9800")
p2 <- plot_distribution(dfsub, "Duration", title = "Total Duration", subtitle = "in minutes", fill = "#F44336")
p3 <- plot_distribution(dfsub, "Break_Duration", title = "Break Duration", subtitle = "in minutes", fill = "#3F51B5")

p4 <- plot_waffle(dfsub, "Sex") +
  scale_fill_manual(values = c("Male" = "#2196F3", "Female" = "#E91E63", "Other" = "#FF9800"))

p5 <- plot_waffle(dfsub, "Education") +
  scale_fill_viridis_d()

p6 <- plot_waffle(dfsub, "Nationality") +
  scale_fill_metro_d()

p7 <- plot_waffle(dfsub, "Ethnicity") +
  scale_fill_manual(values = c("Latino" = "#FF5722", "Asian" = "#FF9800", "Caucasian" = "#2196F3", "African" = "#4CAF50", "Jewish" = "#9C27B0"))

p8 <- plot_waffle(dfsub, "Screen_Resolution") +
  scale_fill_pizza_d()

p9 <- plot_waffle(dfsub, "Device_OS") +
  scale_fill_bluebrown_d()

# p10 <- plot_waffle(dfsub, "Screen_Refresh") +
#   scale_fill_viridis_d()


(p1 / p2 / p3) | (p4 / p5 / p6) | (p7 / p8 / p9)

Results

Delboeuf

Descriptive

plot_descriptive <- function(data, side="leftright") {
  
  if(side == "leftright") {
    x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
    x <- tools::toTitleCase(gsub("arrow", "", x))
    if(x == "Left") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Left", "Right")
    } else if(x == "Right") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Right", "Left")
    }
  } else {
    x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
    x <- tools::toTitleCase(gsub("arrow", "", x))
    if(x == "Up") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Up", "Down")
    } else if(x == "Down") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Down", "Up")
    }
    data$Answer <- fct_rev(data$Answer)
  }
  
  dodge1 <- 0.1 * diff(range(data$Illusion_Difference))
  dodge2 <- -0.1 * diff(range(data$Illusion_Strength))
  
  colors <- colorRampPalette(c("#4CAF50", "#009688", "#00BCD4", "#2196F3", "#3F51B5", "#673AB7", "#9C27B0"))(length(unique(data$Illusion_Strength)))
  
  p1 <- data |> 
    group_by(Illusion_Difference, Illusion_Strength, Answer) |> 
    summarize(Error = mean(Error)) |> 
    mutate(Illusion_Strength = as.factor(round(Illusion_Strength, 2))) |> 
    ggplot(aes(x = Illusion_Difference, y = Error)) +
    geom_bar(aes(fill=Illusion_Strength), position = position_dodge(width=dodge1), stat="identity") +
    geom_line(aes(color = Illusion_Strength), position = position_dodge(width=dodge1)) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_manual(values = colors) +
    scale_fill_manual(values = colors) +
    theme_modern() +
    labs(
      color = "Illusion Strength", 
      fill = "Illusion Strength",
      y = "Probability of Error",
      x = "Task Difficulty"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
  
  colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data$Illusion_Difference)))
    
  p2 <- data |> 
    group_by(Illusion_Difference, Illusion_Strength, Answer) |> 
    summarize(Error = mean(Error)) |> 
    mutate(Illusion_Difference = as.factor(round(Illusion_Difference, 2))) |> 
    ggplot(aes(x = Illusion_Strength, y = Error)) +
    geom_vline(xintercept=0, linetype="dotted", alpha=0.6) +
    geom_bar(aes(fill=Illusion_Difference), position = position_dodge(width=dodge2), stat="identity") +
    geom_line(aes(color = Illusion_Difference), position = position_dodge(width=dodge2)) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_manual(values = colors) +
    scale_fill_manual(values = colors) +
    theme_modern() +
    labs(
      color = "Task Difficulty", 
      fill = "Task Difficulty",
      y = "Probability of Error",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5)) 
  
  if(side == "leftright") {
    p <- ((p1 + facet_wrap(~Answer, ncol=2, labeller = "label_both")) /
      (p2 + facet_wrap(~Answer, ncol=2, labeller = "label_both"))) + 
    plot_annotation(title = paste(data$Illusion_Type[1], "Illusion"), 
                    theme = theme(plot.title = element_text(face = "bold", hjust = 0.5)))
  } else {
    p <- ((p1 + facet_wrap(~Answer, nrow=2, labeller = "label_both")) |
      (p2 + facet_wrap(~Answer, nrow=2, labeller = "label_both"))) + 
    plot_annotation(title = paste(data$Illusion_Type[1], "Illusion"), 
                    theme = theme(plot.title = element_text(face = "bold", hjust = 0.5)))
  }
  p
}

data <- filter(df, Illusion_Type == "Delboeuf")

plot_descriptive(data)

Model Selection

best_models <- function(data) {
  models <- list()
  for(i in 1:1) {
    for(j in 1:1) {
      for(k1 in c("", "_log", "_sqrt", "_cbrt")) { 
        for(k2 in c("")) {  
          for(side in c("", "-side")) {
            name <- paste0("dif", k1, i, "-", "str", k2, j, side)
            # print(name)
            f <- paste0("poly(Illusion_Difference", 
                        k1,
                        ", ",
                        i,
                        ") * poly(Illusion_Strength",
                        k2, 
                        ", ",
                        j, 
                        ") + (1|Participant)")
            
            if(side == "-side") f <- paste0("Illusion_Side * ", f)
            
            m <- glmmTMB::glmmTMB(as.formula(paste0("Error ~ ", f)), 
                                  data = data, family = "binomial")
            if(performance::check_convergence(m)) {
              models[[name]] <- m
            }
          }
        }
      }
    }
  }

  to_keep <- compare_performance(models, metrics = c("BIC")) |> 
    arrange(BIC) |> 
    slice(1:10) |> 
    pull(Name)
  
  
  test <- test_performance(models[to_keep], reference=1)
  perf <- compare_performance(models[to_keep], metrics = c("BIC", "R2")) 
  
  merge(perf, test) |> 
    arrange(BIC) |> 
    select(Name, BIC, R2_marginal, BF) |> 
    mutate(BF = insight::format_bf(BF, name=""))
}

best_models(data)
##                  Name  BIC R2_marginal       BF
## 1      dif_cbrt1-str1 3305       0.392         
## 2      dif_sqrt1-str1 3307       0.400  = 0.353
## 3       dif_log1-str1 3317       0.412  = 0.003
## 4 dif_cbrt1-str1-side 3330       0.394  < 0.001
## 5 dif_sqrt1-str1-side 3332       0.401  < 0.001
## 6           dif1-str1 3333       0.424  < 0.001
## 7  dif_log1-str1-side 3341       0.414  < 0.001
## 8      dif1-str1-side 3356       0.426  < 0.001

Model Visualization

cbrt <- function(x) sign(x) * abs(x)**(1/3)

formula <- brms::bf(
  Error ~ cbrt(Illusion_Difference) * Illusion_Strength + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 12.5 seconds.
## Chain 3 finished in 14.5 seconds.
## Chain 1 finished in 15.4 seconds.
## Chain 2 finished in 16.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 14.6 seconds.
## Total execution time: 16.2 seconds.

# parameters::parameters(model)
plot_model <- function(data, model) {
  data <- mutate(data, .dots_side = ifelse(Error == 1, "bottom", "top"))
  
  # Get variables
  vars <- insight::find_predictors(model)$conditional
  vardiff <- vars[1]
  varstrength <- vars[2]
  
  # Get predicted
  pred <- estimate_relation(model,
                            at = vars,
                            length = c(NA, 25))
  pred[[vardiff]] <- as.factor(pred[[vardiff]])
  
  # Set colors for lines
  colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data[[vardiff]])))
  diffvals <- as.numeric(as.character(unique(sort(pred[[vardiff]]))))
  names(colors) <- diffvals
  
  # Assign color from the same palette to every observation of data (for geom_dots)
  closest <- diffvals[max.col(-abs(outer(data[[vars[1]]], diffvals, "-")))]
  data$color <- colors[as.character(closest)]
  data$color <- fct_reorder(data$color, closest)
  
  # Manual jittering
  xrange <- 0.05*diff(range(data[[varstrength]]))
  data$x <- data[[varstrength]]
  data$x[data$x > 0] <- data$x[data$x > 0] - runif(sum(data$x > 0), 0, xrange)
  data$x[data$x < 0] <- data$x[data$x < 0] + runif(sum(data$x < 0), 0, xrange)
  data$x[round(data$x, 2) == 0] <- data$x[round(data$x, 2) == 0] + runif(sum(round(data$x, 2) == 0), -xrange/2, xrange/2)
  
  
  pred |>
    ggplot(aes_string(x = varstrength, y = "Predicted")) +
    geom_dots(
      data = data,
      aes(x=x, y = Error, group = Error, side = .dots_side, order=color), 
      fill = data$color,
      color = NA, 
      alpha=0.5) +
    geom_slab(data=data, aes(y = Error)) +
    geom_ribbon(aes_string(ymin = "CI_low", ymax = "CI_high", fill = vardiff, group = vardiff), alpha = 0.2) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_hline(yintercept = c(0.05, 0.5, 0.95), linetype = "dotted", alpha = 0.5) +
    geom_line(aes_string(color = vardiff, group = vardiff)) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_manual(values = colors) +
    scale_fill_manual(values = colors) +
    coord_cartesian(xlim=c(min(data[[varstrength]]), max(data[[varstrength]]))) +
    theme_modern() +
    labs(
      title = paste0(data$Illusion_Type[1], " Illusion"),
      color = "Difficulty", fill = "Difficulty",
      y = "Probability of Error",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

plot_model(data, model)

GAM

make_gam <- function(data) {
  
  formula <- brms::bf(
    Error ~ t2(Illusion_Difference, Illusion_Strength, bs = "cr", k=4) + 
      (1 | Participant),
    family = "bernoulli"
  )

  model <- brms::brm(formula,
    data = data,
    refresh = 0
  )
  
  list(p = plot_model(data, model), model = model)
}

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 28.2 seconds.
## Chain 3 finished in 28.5 seconds.
## Chain 1 finished in 28.9 seconds.
## Chain 2 finished in 30.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 28.9 seconds.
## Total execution time: 30.2 seconds.
gam$p

# Parameter Standardization

std_params <- function(model, min=0, max=2) {
  estimate_relation(
  model,
    at = list(Illusion_Strength = c(0), 
              Illusion_Difference = seq(min, max, length.out=500)),
    ) |> 
    select(Illusion_Strength, Illusion_Difference, Error = Predicted) |> 
    slice(c(which.min(abs(Error - 0.005)), 
            which.min(abs(Error - 0.025)), 
            which.min(abs(Error - 0.25)))) |> 
    mutate(Error = insight::format_value(Error, as_percent=TRUE))
}

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=2)
std_params(gam$model, min=0, max=2)

Ebbinghaus

Descriptive

data <- filter(df, Illusion_Type == "Ebbinghaus")

plot_descriptive(data)

Model Selection

best_models(data)
##                  Name  BIC R2_marginal       BF
## 1      dif_cbrt1-str1 2212       0.615         
## 2      dif_sqrt1-str1 2215       0.615  = 0.199
## 3 dif_cbrt1-str1-side 2227       0.623  < 0.001
## 4 dif_sqrt1-str1-side 2230       0.623  < 0.001
## 5       dif_log1-str1 2232       0.617  < 0.001
## 6           dif1-str1 2258       0.619  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ cbrt(Illusion_Difference) * Illusion_Strength + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 18.4 seconds.
## Chain 3 finished in 18.6 seconds.
## Chain 2 finished in 18.8 seconds.
## Chain 1 finished in 19.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 18.7 seconds.
## Total execution time: 19.2 seconds.

# parameters::parameters(model)
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 29.6 seconds.
## Chain 4 finished in 30.3 seconds.
## Chain 1 finished in 32.0 seconds.
## Chain 3 finished in 33.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 31.3 seconds.
## Total execution time: 33.2 seconds.
gam$p

# Parameter Standardization

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=2)
std_params(gam$model, min=0, max=2)

Rod and Frame

Descriptive

data <- filter(df, Illusion_Type == "Rod-Frame")

plot_descriptive(data)


data <- filter(data, abs(Illusion_Strength) < 15)

plot_descriptive(data)

Model Selection

best_models(data)
##                  Name  BIC R2_marginal       BF
## 1      dif_sqrt1-str1 2070       0.464         
## 2       dif_log1-str1 2072       0.458  = 0.371
## 3      dif_cbrt1-str1 2075       0.452  = 0.063
## 4 dif_sqrt1-str1-side 2095       0.478  < 0.001
## 5           dif1-str1 2095       0.505  < 0.001
## 6  dif_log1-str1-side 2097       0.471  < 0.001
## 7 dif_cbrt1-str1-side 2101       0.465  < 0.001
## 8      dif1-str1-side 2119       0.525  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * poly(Illusion_Strength, 2) + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 6.4 seconds.
## Chain 3 finished in 6.6 seconds.
## Chain 4 finished in 6.6 seconds.
## Chain 2 finished in 6.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 6.5 seconds.
## Total execution time: 6.9 seconds.

# parameters::parameters(model)
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 13.8 seconds.
## Chain 2 finished in 16.5 seconds.
## Chain 4 finished in 17.1 seconds.
## Chain 3 finished in 17.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 16.2 seconds.
## Total execution time: 17.5 seconds.
gam$p

# Parameter Standardization

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0.01, max=12)
std_params(gam$model, min=0.01, max=12)

Vertical-Horizontal

Descriptive

data <- filter(df, Illusion_Type == "Vertical-Horizontal")

plot_descriptive(data)


data <- filter(data, abs(Illusion_Strength) < 90)

plot_descriptive(data)

Model Selection

best_models(data)
##                  Name  BIC R2_marginal       BF
## 1      dif_sqrt1-str1 2011       0.661         
## 2      dif_cbrt1-str1 2014       0.657  = 0.234
## 3       dif_log1-str1 2017       0.674  = 0.050
## 4           dif1-str1 2022       0.679  = 0.004
## 5 dif_sqrt1-str1-side 2031       0.665  < 0.001
## 6 dif_cbrt1-str1-side 2034       0.662  < 0.001
## 7  dif_log1-str1-side 2038       0.677  < 0.001
## 8      dif1-str1-side 2043       0.681  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ sqrt(Illusion_Difference) * Illusion_Strength + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 22.4 seconds.
## Chain 1 finished in 22.6 seconds.
## Chain 3 finished in 22.8 seconds.
## Chain 2 finished in 23.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 22.9 seconds.
## Total execution time: 24.0 seconds.

# parameters::parameters(model)
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 22.6 seconds.
## Chain 1 finished in 22.7 seconds.
## Chain 4 finished in 23.1 seconds.
## Chain 2 finished in 23.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 23.0 seconds.
## Total execution time: 23.8 seconds.
gam$p

# Parameter Standardization

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=0.40)
std_params(gam$model, min=0, max=0.40)

Zöllner

Descriptive

data <- filter(df, Illusion_Type == "Zöllner")

plot_descriptive(data)


data <- filter(data, abs(Illusion_Strength) < 45)

plot_descriptive(data)

Model Selection

best_models(data)
##                  Name  BIC R2_marginal       BF
## 1      dif_cbrt1-str1 1041       0.405         
## 2       dif_log1-str1 1043       0.414  = 0.464
## 3      dif_sqrt1-str1 1044       0.423  = 0.260
## 4           dif1-str1 1065       0.491  < 0.001
## 5 dif_cbrt1-str1-side 1066       0.409  < 0.001
## 6  dif_log1-str1-side 1068       0.418  < 0.001
## 7 dif_sqrt1-str1-side 1069       0.427  < 0.001
## 8      dif1-str1-side 1090       0.494  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ cbrt(Illusion_Difference) * Illusion_Strength + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 9.8 seconds.
## Chain 4 finished in 10.1 seconds.
## Chain 2 finished in 10.7 seconds.
## Chain 3 finished in 12.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 10.6 seconds.
## Total execution time: 12.0 seconds.

# parameters::parameters(model)
plot_model(data, model)

Parameter Standardization

# Parameter Standardization

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=5)

White

Descriptive

data <- filter(df, Illusion_Type == "White")

plot_descriptive(data)

Model Selection

best_models(data)
##                  Name  BIC R2_marginal       BF
## 1       dif_log1-str1 2176       0.783         
## 2      dif_cbrt1-str1 2178       0.784  = 0.404
## 3      dif_sqrt1-str1 2181       0.784  = 0.101
## 4  dif_log1-str1-side 2188       0.788  = 0.003
## 5 dif_cbrt1-str1-side 2190       0.789  = 0.001
## 6 dif_sqrt1-str1-side 2192       0.790  < 0.001
## 7           dif1-str1 2195       0.787  < 0.001
## 8      dif1-str1-side 2206       0.793  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ log(Illusion_Difference) * poly(Illusion_Strength, 2) + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 79.3 seconds.
## Chain 2 finished in 89.6 seconds.
## Chain 3 finished in 90.0 seconds.
## Chain 4 finished in 93.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 88.0 seconds.
## Total execution time: 93.5 seconds.

# parameters::parameters(model)
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 39.7 seconds.
## Chain 1 finished in 43.7 seconds.
## Chain 4 finished in 51.6 seconds.
## Chain 2 finished in 52.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 46.8 seconds.
## Total execution time: 52.2 seconds.
gam$p

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
# unique(data$Illusion_Strength)

std_params(model, min=0, max=20)
std_params(gam$model, min=0, max=20)

Müller-Lyer

Descriptive

data <- filter(df, Illusion_Type == "Müller-Lyer")

plot_descriptive(data, side = "updown")

Model Selection

best_models(data)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
##                  Name  BIC R2_marginal       BF
## 1      dif_cbrt1-str1 2449       0.749         
## 2  dif_log1-str1-side 2534       0.748  < 0.001
## 3 dif_cbrt1-str1-side 2535       0.757  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ log(Illusion_Difference) * poly(Illusion_Strength, 2) + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 16.7 seconds.
## Chain 4 finished in 18.4 seconds.
## Chain 3 finished in 19.4 seconds.
## Chain 1 finished in 19.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 18.5 seconds.
## Total execution time: 19.7 seconds.

# parameters::parameters(model)
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 32.9 seconds.
## Chain 2 finished in 33.8 seconds.
## Chain 3 finished in 33.8 seconds.
## Chain 1 finished in 34.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 33.6 seconds.
## Total execution time: 34.1 seconds.
gam$p

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=0.6)
std_params(gam$model, min=0, max=0.6)

Ponzo

Descriptive

data <- filter(df, Illusion_Type == "Ponzo")

plot_descriptive(data, side = "updown")

Model Selection

best_models(data)
##                 Name  BIC R2_marginal       BF
## 1     dif_cbrt1-str1 2551       0.661         
## 2     dif_sqrt1-str1 2555       0.665  = 0.187
## 3      dif_log1-str1 2576       0.678  < 0.001
## 4          dif1-str1 2590       0.684  < 0.001
## 5 dif_log1-str1-side 2592       0.687  < 0.001
## 6     dif1-str1-side 2606       0.694  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ poly(Illusion_Difference, 2) * poly(Illusion_Strength, 2) + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 12.3 seconds.
## Chain 2 finished in 13.2 seconds.
## Chain 1 finished in 13.3 seconds.
## Chain 3 finished in 14.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 13.4 seconds.
## Total execution time: 14.9 seconds.

# parameters::parameters(model)
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 29.3 seconds.
## Chain 2 finished in 30.4 seconds.
## Chain 3 finished in 30.3 seconds.
## Chain 4 finished in 35.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 31.3 seconds.
## Total execution time: 35.3 seconds.
gam$p

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=20)
std_params(gam$model, min=0, max=20)

Poggendorff

Descriptive

data <- filter(df, Illusion_Type == "Poggendorff")

plot_descriptive(data, side = "updown")


data <- filter(data, abs(Illusion_Strength) < 45)

plot_descriptive(data, side = "updown")

Model Selection

best_models(data)
##                  Name  BIC R2_marginal       BF
## 1      dif_cbrt1-str1 1900       0.527         
## 2      dif_sqrt1-str1 1905       0.526  = 0.086
## 3       dif_log1-str1 1923       0.522  < 0.001
## 4 dif_cbrt1-str1-side 1929       0.529  < 0.001
## 5           dif1-str1 1930       0.521  < 0.001
## 6 dif_sqrt1-str1-side 1934       0.529  < 0.001
## 7  dif_log1-str1-side 1952       0.527  < 0.001
## 8      dif1-str1-side 1958       0.526  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ cbrt(Illusion_Difference) * poly(Illusion_Strength, 2) + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 11.2 seconds.
## Chain 1 finished in 11.6 seconds.
## Chain 3 finished in 11.8 seconds.
## Chain 4 finished in 11.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 11.6 seconds.
## Total execution time: 12.1 seconds.

# parameters::parameters(model)
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 18.1 seconds.
## Chain 4 finished in 19.0 seconds.
## Chain 3 finished in 19.4 seconds.
## Chain 2 finished in 21.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 19.4 seconds.
## Total execution time: 21.1 seconds.
gam$p

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=0.5)
std_params(gam$model, min=0, max=0.5)

Contrast

Descriptive

data <- filter(df, Illusion_Type == "Contrast")

plot_descriptive(data, side = "updown")

Model Selection

best_models(data)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
##                  Name  BIC R2_marginal       BF
## 1      dif_sqrt1-str1 2353       0.673         
## 2      dif_cbrt1-str1 2354       0.673  = 0.713
## 3       dif_log1-str1 2357       0.672  = 0.128
## 4           dif1-str1 2359       0.675  = 0.047
## 5 dif_cbrt1-str1-side 2580       0.685  < 0.001
## 6      dif1-str1-side 2585       0.685  < 0.001

Model Visualization

formula <- brms::bf(
  Error ~ Illusion_Difference * poly(Illusion_Strength, 2) + 
    (1 | Participant),
  family = "bernoulli"
)

model <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 22.7 seconds.
## Chain 3 finished in 23.0 seconds.
## Chain 1 finished in 23.3 seconds.
## Chain 2 finished in 24.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 23.4 seconds.
## Total execution time: 24.5 seconds.
plot_model(data, model)

GAM

gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 30.7 seconds.
## Chain 1 finished in 31.8 seconds.
## Chain 3 finished in 33.8 seconds.
## Chain 4 finished in 34.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 32.7 seconds.
## Total execution time: 34.6 seconds.
gam$p

# Parameter Standardization

# range(data$Illusion_Difference)
# range(data$Illusion_Strength)

std_params(model, min=0, max=25)
std_params(gam$model, min=0, max=25)